Electronic Quantum Monte Carlo Calculations of Atomic Forces, Vibrations, and 

Anharmonicities 
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Atomic forces are calculated for first-row monohydrides and carbon monoxide within electronic 
quantum Monte Carlo (QMC). Accurate and efficient forces are achieved by using an improved 
method for moving variational parameters in variational QMC. Newton's method with singular 
value decomposition (SVD) is combined with steepest descent (SD) updates along directions re- 
jected by the SVD, after initial SD steps. Dissociation energies in variational and diffusion QMC 
agree well with experiment. The atomic forces agree quantitatively with potential energy surfaces, 
demonstrating the accuracy of this force procedure. The harmonic vibrational frequencies and 
anharmonicity constants, derived from the QMC energies and atomic forces, also agree well with 
experimental values. 



I. INTRODUCTION 

Quantum Monte Carlo (QMC) is an effective method 
for solving the time-independent Schrodinger equation, 
and has become quite successful in computing ground- 
state total energies. The QMC method gives energies of 
atoms, molecules, and solids that are comparably accu- 
rate or more accurate than traditional techniques such 
as density functional theory (DFT), multiconfiguration 
self-consistent field (MCSCF), or coupled cluster meth- 
ods. Although the situation for the calculation of prop- 
erties other than energies has been less favorable, the 
accurate QMC calculation of atomic forces has been en- 
abled through the recent developments made in this area 
by Assaraf and Caffarel Q, , Filippi and Umrigar Q , 
Casalegno, Mella, and Rappe 4], Chiesa, Ceperley, and 
Zhang 5], and others. 

In this paper, we extend our atomic force methodology 
to all the first-row monohydrides and carbon monoxide. 
In order to acquire energies and forces efficiently for these 
systems, we also describe an improved algorithm for op- 
timizing variational Monte Carlo (VMC) wave functions. 
As in our previous paper the first and second deriva- 
tives of the variational energy are analytically computed, 
and used to perform Newton's method parameter up- 
dates with SVD. We now propose augmenting this ap- 
proach by using the steepest descent (SD) method in 
the subspace neglected by the Newton's method with 
SVD. In the initial stage of parameter update, Newton's 
method might give poor result since the second deriva- 
tives include larger noise when the parameters are far 
from the optimum. So we take two SD steps before start- 
ing Newton's method. The improved algorithm was ap- 
plied to the calculation of the ground-state energies and 
forces of the first-row monohydrides and carbon monox- 
ide. In general, the direct application of the variational 
principle yields significantly lower energy than variance 
minimization methods, so minimizing the energy is ad- 
vantageous. The wave functions optimized in VMC were 
used as a guiding function to compute more accurate en- 



ergies and forces in diffusion Monte Carlo (DMC). 

In this paper, total energies, dissociation energies, 
forces, harmonic vibrational frequencies, and anhar- 
monicity constants are reported for all first-row mono- 
hydrides from LiH to HF, as well as for CO. In all cases, 
the computed results agree well with experiment. The 
dissociation energies in VMC are significantly improved 
with respect to a previous VMC study of the hydrides. 



II. THEORETICAL BACKGROUND AND 
COMPUTATIONAL DETAILS 

The variational parameters used in VMC will be de- 
noted as ci, C2, ■ • ■, c„, here. The VMC energy expec- 
tation value, Et, is a function of these variational pa- 
rameters, and the parameter set that minimizes Ej^ is 
sought. 

The SD method is useful in the initial stages of param- 
eter optimization in VMC, due to the large error bars of 
Hessian matrix components. One arbitrary constant is 
necessary to implement the SD method. We used the 
following two-step scheme to find a good SD constant. 

Let Qo and Q be the vectors composed of variational 
parameters before and after update, respectively: 



Qo = (ci,o c 2 , 
Q = (ci c 2 • • 



Cn,0j 



(1) 

(2) 



And let g be the gradient vector of energy with respect 
to the variational parameters: 

m ••• 9 n) T =m n ••• m T - (3) 

In the first update, a value, a^°\ is chosen as a SD con- 
stant, which is small enough not to exhaust the downhill 
direction, 



Q = Qo a (0) g(Qo) 



(4) 
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After the first update, VMC simulation is performed 
again to get the gradient at the new parameter set, Q. 
If we consider only the ith component, the best value for 
the new SD constant, a\ , will make the gradient com- 
ponent, <7i, zero in the next simulation, and will be given 
by 



,(0) 



i-<7i(Q)/<?i(Qo)' 



(5) 



Although values are different from component to 
component, it is usually the case that they are quite sim- 
ilar. So the averaged value was used for the next update: 



1 ™ 

n ^ 

i=i 



,(!) 



(G) 



These two steps of parameter updates using SD reduce 
the energy enough to greatly reduce the error bars, en- 
abling the Newton's method. 

If we let H represent the Hessian matrix, the parame- 
ters can be updated according to Newton's method, 



Q = Qo-H- 1 (Q )g(Qo). 



(7) 



Since H(Q ) and g(Qo) are calculated in the VMC sim- 
ulation, we must invert H(Qo) for the Newton's method 
update of the parameters. 

It is well-known that any matrix, e.g., H, can be ex- 
pressed as 



H = U[diag( Wj )]V T , 



(8) 



where Wj > and U,V are orthogonal Q. For a square 
matrix, the inverse matrix can be obtained by 



H" 1 = V[diag(l/^)]U T . 



(9) 



Very small values of utj lead to erroneous moves along 
the directions corresponding to these components due to 
large 1/wj terms. For that reason, if Wj is less than a 
certain threshold value, 1/uij is set to in the actual 
calculation (SVD). 

SVD method has been tested for the inversion of Hes- 
sian matrix and it consistently gives robust results in 
many cases. However, the SVD method, by zeroing out 
small Wj values, is equivalent to abandoning the corre- 
sponding search directions, the use of which may give 
better result. So we propose a modified algorithm in 
which the SD method is added for components discarded 
in the SVD method. 

If we let U and V be equal to the square matrix whose 
column vectors are the normalized eigenvectors of H, 
{wj} will be the eigenvalues of H. For values of Wj 
that are smaller than the threshold, 1/wj can be replaced 
by a constant, a, instead of zero, which is equivalent to 
the SD method along the corresponding directions. This 
method makes it possible to use the information for all 



directions, some of which are discarded in SVD method, 
and it can be beneficial in cases where some eigenvalues 
of the Hessian matrix become close to zero, due to the 
noise inherent in QMC. In case of SVD algorithm, Wj 
is always nonnegative, which corresponds to the abso- 
lute value of eigenvalue of H. If any eigenvalue is nega- 
tive and its absolute value is larger than the threshold, 
there is a problem that the direction corresponding to 
this is not discarded, even though this does not happen 
so frequently. This small problem of negative eigenvalues 
can be handled by using the modified method with the 
same positive threshold and zero steepest descent con- 
stant, and we used this modified method in the actual 
implementation. 

To construct the trial wave functions used in VMC, the 
following method was used. First, a contracted Gaussian- 
type function (CGTF) was fitted to each Slater-type or- 
bital (STO). Ten primitive Gaussians were used for Is, 
eight for 2s or 2p, and six for 3s, 3p, or 3d type STOs. 
The orbital exponents of STOs in the works of Cade and 
Huo Hi were adopted (excluding the /-type orbitals). 
In case of the first-row monohydrides, each first-row atom 
has 29 STOs centered on it (Is, Is', 2s, 2s', 3s, three 2p's, 
three 2p"s, three 2p"'s, three 3p's, six 3(fs, and six 3d"s 
for Li, and Is, Is', 2s, 2s', 3s, three 2p's, three 2p"s, 
three 2p"'s, three 2p""s, six 3efs, and six 3d"s for other 
first-row atoms) and hydrogen atom has 6 STOs centered 
on it (Is, Is', 2s, and three 2p , s) as a basis set. In case 
of carbon monoxide, each atom has 19 STOs centered on 
it (Is, Is', 2s, 3s, three 2p , s, three 2p"s, three 2p'"s, and 
six 3g?'s) as a basis set. 

Each molecular orbital (MO) was expressed as a lin- 
ear combination of STOs, the coefficients of which were 
obtained using the Hartree-Fock method in Gaussian 98 
(G98) [lCj. For the open shell molecules, restricted open 
shell Hartree-Fock (ROHF) wave functions were used. 
The MOs from G98 were used to construct the Slater 
determinants for a and (3 electrons. While multidetermi- 
nant trial wave function gives improved results for some 
systems, it was reported that the use of single determi- 
nant trial wave function gave good results in the calcu- 
lations of the first-row hydrides [HIH. Since the use 
of multideterminant trial wave function is much more 
time-consuming, we used only single determinant in the 
calculation here. The product of two determinants was 
multiplied by a positive correlation factor to form a trial 
wave function 



a positr 



a,i<i 



where 



Uaij — c ka{Tai 



aj ai J ij 



(10) 



(11) 



In this equation, a and i,j refer to the nuclei and 
the electrons, respectively, and f is defined by f = 
br/(l + br). Cfea's are variational Jastrow parameters. 



3 



We used b 



1 a n and included 30 terms for di- 



atomic molecules, namely, 4 electron-electron, 6 electron- 
nucleus, and 20 electron-electron-nucleus terms. In case 
of atoms, we used 17 parameters composed of 4 electron- 
electron, 3 electron- nucleus, and 10 electron-electron- 
nuclcus terms, to be consistent with the calculation of 
diatomic molecules. 

Five different bond distances around the experimental 
bond length were used for calculation, namely 90%, 95%, 
100%, 105% and 110% of the experimental bond length, 
r exp . 2000 walkers were used for all the calculations in 
this paper. In updating Jastrow parameters, average over 
100 blocks was made typically, where each block was the 
average over 100 steps. To accelerate the sampling, a 
Fokker-Planck type equation was used 

After a short initial simulation without Jastrow fac- 
tor, the Hartree-Fock wave function was multiplied by 
the Jastrow factor with all parameters set to zero. The 
gradient and Hessian of energy with respect to the Jas- 
trow parameters were computed in the VMC simulation 
after this step. Using the gradient and Hessian informa- 
tion, a new Jastrow parameter set is calculated, and a 
new VMC simulation is performed with this updated pa- 
rameter set. This process was iterated until the energy 
converged. Fully optimized parameters were obtained by 
10-15 iterations. One iteration took about 30 minutes for 
LiH and about 90 minutes for HF when a single 2.8 GHz 
Intel® Xeon Processor was used. 

After optimizing the trial wave function using VMC, 
a fixed-node DMC calculation was performed using im- 
portance sampling, as proposed by Reynolds, Ceperley, 
Alder, and Lester The DMC time step was 0.005 

a.u. for the first-row hydrides and 0.0005-0.001 a.u. for 
carbon monoxide. A similar DMC method was used by 
Liichow and Anderson [Tll[l7T | in their calculation of first- 
row hydrides. 

Force calculations were performed in both VMC and 
DMC. We followed the method described previously £|. 
If the wave function were exact, the exact force would be 
given by the Hellmann-Feynman theorem (HFT). Since 
the trial wave function, is not exact, terms that can- 
cel in case of exact wave functions should be considered, 
in addition to the HFT expression. Retaining terms in- 
volving wave function derivatives gives the total atomic 
force on atom a in direction q: 



F„, 



pHFT 

qa 



pPulay , pc 
qa r x ga' 



where 



pPulay 
qa 

and 



F, 



\c)R a , 



HFT 



(*7 



I dH 
1 dR„a 



1*7 



(*T *t) 



2(E) 



\ dR„, 



1*7 



VMC" 



;*t|*t! 



(12) 



(13) 



-, (14) 



F; 



qa 



E 



dc k d(E) YMC 
dR qa dc k 



(15) 



These expressions apply for VMC, and similar equations 
are used for DMC simulations |4|- F^ lay incorporates 
the explicit dependence of the wave function on the nu- 
clear coordinates (Pulay's correction ^H), an< 4 can ^ e 
easily calculated through VMC or DMC simulations. F^ a 
depends implicitly on the nuclear coordinates through 
the variational parameters. However, since an energy- 
minimized wave function is used, i.e., d(E)yMc/dck = 0, 
this force term makes zero contribution. In the calcula- 
tion of the Hellmann-Feynman theorem force, F^ FT , the 
rcnormalizcd estimator proposed by Assaraf and Caffarel 
was used to reduce the variance of the force calcu- 
lation. The expectation value of this estimator, F^p , 
is the same as F^ FT , but the variance of the former is 
much smaller. In our force calculation, F^ + F^ lay 
was computed by averaging over the walkers. 



III. RESULTS AND DISCUSSION 

The energies of first-row monohydrides and carbon 
monoxide at various bond distances were calculated. The 
plot of energy versus bond distance for hydrogen fluoride 
(HF) is shown in Figure^ In obtaining each point, 1000 
blocks, each of which was composed of 100 steps, were 
used with optimized Jastrow parameters. The plots for 
other molecules are similar to that for HF. The energies 
obtained from VMC are a few tenths of a Hartree lower 
than the Hartree-Fock energies obtained from G98, so 
the Hartree-Fock results are not shown in the figure. It 
can be seen from Table [|] that the DMC energy is signif- 
icantly lower than the VMC energy and is close to the 
experimental value. 

The bond dissociation energies, D e , were calculated by 
taking the differences between QMC energies of diatomic 
molecules in Table Q] and QMC energies of atoms. To 
be consistent in the number of Jastrow parameters, we 
did the calculation of atoms with 17 parameters. The 
VMC energies of atoms with 17 parameters falls between 




1.6 1.7 1.8 1.9 

Bond Distance (Bohr) 

FIG. 1: Energy and force calculation of HF with VMC and 
DMC. Two thin horizontal lines at each data point show the 
energy error bar. The slope of the thick lines show the force 
at each data point. 
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TABLE I: Eq, r e , u e , and uj e x e for LiH - HF and CO obtained 
from VMC and DMC calculation and experimental data. 







E (Ha) 


r e (Bohr) 


u e (cm -1 ) 


LU e X e (cm 1 ) 


LiH 


VMC 


-8.063 


3.038(1) 


1402(4) 


25.7(1) 




DMC 


-8.070 


3.020(1) 


1417(4) 


24.8(1) 




Exp 


-8.070 


3.015 


1406 


23.2 


BeH 


VMC 


-15.235 


2.519(1) 


2141(4) 


56.6(2) 




DMC 


-15.246 


2.515(1) 


2134(4) 


58.5(2) 




Exp 


-15.248 


2.537 


2061 


36.3 


BH 


VMC 


-25.254 


2.370(1) 


2332(5) 


47.0(2) 




DMC 


-25.275 


2.386(1) 


2369(5) 


47.3(2) 




Exp 


-25.289 


2.329 


2367 


49.4 


CH 


VMC 


-38.438 


2.097(1) 


2961(6) 


77.2(3) 




DMC 


-38.463 


2.112(1) 


2898(6) 


71.8(3) 




Exp 


-38.490 


2.116 


2858 


63.0 


NH 


VMC 


-55.178 


1.941(1) 


3415(7) 


104.3(4) 




DMC 


-55.206 


1.962(1) 


3253(7) 


92.0(4) 




Exp 


-55.247 


1.958 


3282 


78.4 


OH 


VMC 


-75.687 


1.820(1) 


3854(7) 


101.2(4) 




DMC 


-75.720 


1.843(1) 


3690(7) 


91.4(4) 




Exp 


-75.778 


1.832 


3738 


84.9 


HF 


VMC 


-100.407 


1.729(1) 


4206(9) 


89.9(4) 




DMC 


-100.442 


1.755(1) 


4040(9) 


82.4(4) 




Exp 


-100.531 


1.733 


4138 


89.9 


CO 


VMC 


-113.176 


2.095(1) 


2539(16) 


21.1(3) 




DMC 


-113.286 


2.116(2) 


2251(26) 


14.2(3) 




Exp 


-113.377 


2.132 


2170 


13.3 



those with 9 parameters and those with 42 parameters 
reported in Ref. [||. The dissociation energies are sum- 
marized in Tabled together with the results given in the 
work by Liichow and Anderson . Our VMC dissocia- 
tion energies are much closer to experimental values than 
those given by Liichow and Anderson, while our DMC 
results are quite similar to theirs. The improvement in 
our VMC result may be attributed to the effectiveness 
of energy minimization method relative to the variance 
minimization method used for VMC calculations in Ref. 

while part of the improvement is also due to the 
larger number of Jastrow parameters in our calculation. 

Energies calculated by DMC are quite close to the 
experimental values for lighter first-row hydrides, while 
slightly higher energies than experimental values are ob- 
tained for heavier molecules. This may be due to the 
approximations used in DMC calculations: fixed node 
approximation, neglect of the relativistic effect, and the 
error related with finite time step. To estimate the finite 
time step error, DMC calculations at r cxp with several 
different time step values ranging from 0.0001 to 0.005 
a.u. were carried out for first-row hydrides. All calcu- 
lated energies agreed within 2-3 mHartree. 



TABLE II: Dissociation energies D e in kcal/mol for LiH - 
HF and CO obtained from our QMC calculation and from 
literature. 





VMC a 


DMC 


VMC 6 


DMC 6 


Exp 6 


LiH 


54.7 


57.8 


45.7 


57.8 


58.0 


BeH 


57.9 


55.7 


49.4 


52.1 


49.8 


BH 


82.7 


84.7 


63 


84.8 


84.1 


CH 


81.1 


83.5 


81 


83.9 


83.9 


NH 


80.2 


82.3 


77 


81.4 


80.5-84.7 


OH 


105.1 


106.4 


86 


106.4 


106.6 


HF 


140.4 


141.4 


130 


141.3 


141.5 


CO 


218.1 


254.9 






258.7 



"Differences between QMC energies of molecules in Table H] and 
QMC energies of atoms calculated with 17 parameters. 

''From Ref. 1111 for first-row hydrides. 

In the VMC calculation of HF, the Jastrow parameter 
set at r eX p was optimized first, and after the optimization 
at this distance, the bond distance was changed, and the 
MO coefficients corresponding to this bond distance were 
introduced. Then, the Jastrow parameters were reopti- 
mized at this new bond distance. This method makes 
it possible to reduce the CPU time for the calculation at 
other bond distances once the parameter set is optimized 
at one bond distance. This approach is effective because 
the Jastrow parameter sets at different bond distances 
can be quite similar, as measured by the cosine similar- 
ity |19| between Jastrow parameter sets, 



vQm ' Qm\/Qn ' On 

which is close to unity if two vectors are similar. This 
is certainly the case for Jastrow parameter sets of HF 
at various bond distances, as shown in Table ITTT1 This 
approach seems to be useful for the molecular dynamics 
(MD) simulation coupled with QMC, proposed by Gross- 
man and Mitas |2Cj. On the other hand, in case of CH, 
NH, or OH, it was problematic to apply this method and 
we had to optimize the parameters from the beginning 
for all bond distances. The cosine similarity values in 
case of CH are shown in Table ITTA when the parameters 
are optimized separately from scratch for all bond dis- 
tances. If the parameters of HF at each bond distance 
are optimized from scratch, the cosine similarity values 
are around 0.9 for parameter sets optimized at different 
bond distances, and similar energies can be obtained with 
different sets of parameters. 

The energy of BH at r cxp at various stages of param- 
eter optimization is shown in Figure [21 If the SD steps 
are used for initial stages of parameter optimization (B), 
Newton's method with SVD converges to the lowest en- 
ergy after several iterations. If the initial SD steps are 
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TABLE III: The cosine similarity values between Jastrow pa- 
rameter sets obtained from VMC calculations of HF and CH. 



cos^ (HF) 


0.90 r cxp 


0.95 r cxp 


1.00 r cxp 


1.05 r cxp 


1.10 r cxp 


0.90 r cxp 


1.000 










0.95 r cxp 


0.998 


1.000 








1.00 r cxp 


0.997 


1.000 


1.000 






1.05 r cxp 


0.998 


0.997 


0.997 


1.000 




1.10 Texp 


0.997 


0.999 


0.999 


0.997 


1.000 


cos 6 (CH) 


0.90 r cxp 


0.95 r cxp 


1.00 r cxp 


1.05 r cxp 


1.10 Texp 


0.90 r cxp 


1.000 










0.95 r cxp 


0.836 


1.000 








1.00 r cxp 


0.842 


0.964 


1.000 






1.05 r cxp 


0.658 


0.807 


0.701 


1.000 




1.10 r cxp 


0.829 


0.879 


0.848 


0.817 


1.000 



not used (A), Newton's method is somewhat difficult to 
apply due to the large error bars of Hessian components. 
In this case, it was necessary to set the SVD threshold 
somewhat high and to calculate for a long period of time. 
Within this approach, using only the Newton's method 
with SVD does not yield fully optimized energy. The 
simultaneous application of Newton's method and SD 
(steps 6-9) was very useful in this case for more thor- 
ough minimization. 

Forces were computed for each monohydride and car- 
bon monoxide at each bond length in VMC and DMC. 
The force result for HF is shown in Figure where the 
slopes of the line segments superimposed on the energy 
result correspond to the negative of the calculated forces. 
The calculated forces of HF are shown in Table IIVI to- 
gether with the values obtained by fitting energy result to 
the parabolic potential and then calculating the slopes. 
The force at 0.90 r exp is larger than the magnitude of the 
slope of the parabola obtained from the energy result, 



-25.10 



-25.15 



a 
X 



a 

m 



■25.20 



-25.25 



-25.30 




4 6 
Iteration Number 



FIG. 2: The energy of BH at r exp at various stages of parame- 
ter optimization. (A) Newton's method for 1-6 and Newton's 
method with SD for 6-9. (B) Initial SD for 1-3 and Newton's 
method for 3-9. 



while the magnitude of the force at 1.10 r eX p is smaller 
than the parabola tangent, which clearly shows the de- 
viation of the calculated forces from harmonic behavior 
due to anharmonicity. 

The approximate shape of the anharmonic potential 
can be described by the Morse potential pH ]. 



V(r) = D e 1 



,-0(r-r B ) 



(17) 



and this was used in the fitting of the QMC results to 
calculate the properties of diatomic molecules. The vih 
energy level of the Morse potential with reduced mass /i 
is 



E 
he 











(v-\ 


-I) 


- LU e X e I V H 













(18) 



where the harmonic vibrational frequency is given by 
uj e = (3 (100 D e h/2Tr 2 c/i) 1 ' ' 2 and the anharmonicity con- 
stant by uj e x e — (100 h (i 2 /8ir 2 fic) . In this equation, u> e , 
D e and {3 have the unit of cm -1 and other constants are 
in SI units. 

Since we performed QMC calculations at small num- 
ber of bond distances, it would be advantageous to use 
the energy and force results simultaneously in the fitting 
to the Morse potential, which can be accomplished by 
minimizing the following merit function: 



N 



i=l 



Ek-E(n',eL) 



N 

E 



F t -F(r i; a) 

(TFA 



(19) 

Here a is a parameter vector whose components are D e , 
(3, and r e , and the following functional forms were used: 



£7(r;a) = D e 



fl _ e -/3(r-r.)V 



- 1 



(E A + E B ), 

(20) 

F(r; a) = -2D e (3 (l - e-^-*)) e -P^ c ) ^ 

The expression for the energy has been modified to pro- 
duce correct dissociation energy, D e , and Ea and Eb are 
VMC or DMC energies of atoms A and B. The equilib- 
rium bond lengths (r e ), harmonic vibrational frequencies 
(u) e ), and anharmonicity constants (uj e x e ) for all first-row 
monohydride molecules and carbon monoxide calculated 
by fitting energy and force results are summarized in Ta- 
ble HI along with the experimental data The experi- 
mental energies are corrected by adding zero point ener- 
gies. Our calculations agree well with the experimental 
results. 

Each energy or force data point has an error bar as- 
sociated with it, so we followed a simple procedure to 
estimate how the calculated error bars translate into un- 
certainty in other quantities such as equilibrium bond 
length, harmonic vibrational frequency, and anharmonic- 
ity constant. A large set of synthetic data points were 
stochastically generated, such that the average value at 
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TABLE IV: Forces obtained from the slope of parabolic potential energy fits and from the direct calculation for HF. 



Force 


0.90 r oxp 


0.95 r cxp 


1.00 r cxp 


1.05 7-exp 


1.10 r cxp 


VMC (parabola) 
VMC (direct) 


0.113(11) 
0.147(1) 


0.057(6) 
0.060(1) 


0.001(3) 
-0.002(1) 


-0.055(6) 
-0.050(1) 


-0.111(11) 
-0.076(2) 


DMC (parabola) 
DMC (direct) 


0.110(4) 
0.168(1) 


0.056(2) 
0.077(1) 


0.002(1) 
0.015(1) 


-0.051(2) 
-0.033(1) 


-0.105(4) 
-0.064(1) 



each bond length agrees with that obtained from QMC 
with the standard deviation the same as the error bar 
given by the QMC calculation. By computing the av- 
erages and standard deviations of the equilibrium bond 
length (r e ), harmonic vibrational frequency (ui e ), and an- 
harmonicity constant (u) e Xe) for the synthetic data sets, 
the error bars of r e , u> e and uj e x e could be estimated. The 
error bars of the last digit thus calculated are shown in 
parentheses. 



addition of steepest descents to the initial steps and to 
the subspace neglected by Newton's method with SVD 
seems to be advantageous for the molecular systems we 
investigated. 

We could calculate accurate harmonic vibrational 
frequencies and anharmonicity constants of diatomic 
molecules by fitting QMC results to the Morse potential, 
achieving excellent agreement between QMC calculations 
and experiment for these vibrational parameters. 



IV. CONCLUSIONS 

The force calculation method combining energy min- 
imization, Pulay's corrections, and a renormalized 
Hellmann-Feynman estimator worked well with all the 
first-row hydride molecules and carbon monoxide with 
small extra effort. 

The energy minimization method in VMC is useful, 
but it requires an effective optimization scheme. The 
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